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Abstract 

Structural behavior and equation of state of atomic and molecular crystal phases of dense hy- 
drogen at pressures up to 3.5 TPa are systematically investigated with density functional theory. 
The results indicate that the Vinet EOS model that fitted to low-pressure experimental data over- 
estimates the compressibility of dense hydrogen drastically when beyond 500 GPa. Metastable 
multi-atomic molecular phases with weak covalent bonds are observed. When compressed beyond 
about 2.8 TPa, these exotic low-coordinated phases become competitive with the groundstate and 
other high-symmetry atomic phases. Using nudged elastic band method, the transition path and 
the associated energy barrier between these high-pressure phases are evaluated. In particular for 
the case of dissociation of diatomic molecular phase into the atomic metallic Cs-IV phase, the exis- 
tent barrier might raise the transition pressure about 200 GPa at low temperatures. Plenty of flat 
and broad basins on the energy surface of dense hydrogen have been discovered, which should take 
a major responsibility for the highly anharmonic zero point vibrations of the lattice, as well as the 
quantum structure fluctuations in some extreme cases. At zero pressure, our analysis demonstrates 
that all of these atomic phases of dense hydrogen known so far are unquenchable. 
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I. INTRODUCTION 



Having only a single electron outside the nucleus, hydrogen is the simplest and most 
abundant element in the universe. It is also an essential element for models of stellar and 
planetary interiors.- 1 ^ Hydrogen shows characteristics of both the group I alkalis and the the 
group VII halogens. At low pressures, hydrogen isotopes are halogenous, covalent diatomic 
molecules that form insulators. Yet at high pressures, it is one of the most difficult to under- 
stand. It displays anomalous melting behavior with a maximum in the melting temperature 
versus pressure curve at high temperatures,- 1 ^ and undergoes a first-order liquid-liquid tran- 
sition under further compression.-^ At low temperatures, it is experimentally known that 
hydrogen can exist as a rotational crystal (phase I) on a hexagonal close-packed (HCP) 
lattice to high pressures (P < 110 GPa), followed by a transition into the broken-symmetry 
phase II (110 GPa < P < 150 GPa) which is marked by a change to wide-angle libration 
and hence to a continuing incoherence of motion between different molecules, and then to 
phase III at about 150 GPa. 7 Possible pressure-induced insulator-metal transition also has 
been extensively studied up to 320 GPa.- - — However, beyond the fact that protons remain 
paired within this pressure range, their time-average locations are to date experimentally un- 
known, mainly because hydrogen atoms scatter X-rays only weakly, leading to low-resolution 
diffraction patterns.-*^ Experimental data at higher pressures are scarce, and insightful un- 
derstanding of structural behavior of ultra-dense hydrogen is lacking. 

Theoretical prediction of stable crystalline structures and properties of dense hydrogen 
at high pressures has been pursued for decades.—"— It is difficult because of the need to 
search the very large space of possible structures, and the necessity of obtaining accurate 
energies for each of these structures.— 1 ^ First-principles density functional theory (DFT) 
has been proved as an efficient approach of calculating quite accurate energies, and has 
provided insights into properties of various materials, including solid hydrogen under com- 
pressions. At present, DFT offers a high level of theoretical description at which we can 
carry out searches over many possible candidate structures.—"— Recent DFT calculations 
have predicted that within the static-lattice approximation, the most stable phases of dense 
hydrogen are P6 3 /m (<105GPa), C2/c (105-270 GPa), Cmca-12 (270-385 GPa), and Cmca 
(385-490 GPa), followed by atomic I^jamd (Cs-IV) phase,^^ and then R3m (or ffim)& 

On the other hand, at pressures high enough so that electrons are fully ionized out to form 
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a uniform background, hydrogen becomes one-component plasma.— In this regime the domi- 
nant Coulomb interaction forces the ions into an ordered configuration called Wigner crystal, 
and stabilizes in a body-centered cubic (BCC) structure.—"— Except this, the general struc- 
tural and compressional behavior of dense hydrogen at ultra-high pressures beyond several 
Mbar (1 Mbar=100 GPa) are still poorly understood, in spite of a few theoretical investiga- 
tions available.— 1 ^ In particular, the energy barriers of pressure-induced phase transitions 
are completely unknown, especially that of molecule dissociation. A high energy barrier 
would lead to a hysteresis and raise the transition pressure from what it otherwise would 
be.— Another important issue that has not been solved completely is about the applicability 
of extrapolating the equation of state (EOS) that fitted at low pressure to higher pressures. 
This is an intriguing problem because we knew that the EOS model fitted to lower-pressure 
data (up to 42 GPa) 7 fails to capture (underestimates) the compressibility of H2 and D2 at 
the relatively higher pressures,— which is due to the transition of H 2 (or D 2 ) from a freely 
rotating phase into a broken-symmetry phase. At ultra-high pressures, diatomic hydrogen 
dissociates into atomic phase, and this transition might invalidate the early established EOS 
model again very likely. 

The groundstate structures of atomic metallic hydrogen from 500 GPa to 4.5 TPa has 
been extensively searched by McMahon et ai— using the ab initio random structure search- 
ing (AIRSS) method.— Their work was carried out for unit cells containing only 4 and 6 
atoms. However, as recent investigations suggested, complex structures can be adopted by 
simple alkali metals at high pressures.— >22. Similar phenomenon might also occur in dense 
hydrogen. Furthermore, there is no absolute guarantee for a specific structure searching 
approach to find the true ground state within a finite computational time due to the limited 
phase space it can explore. There are examples that AIRSS calculations failed to detect 
lower energy states which later were captured by other structure searching method such as 
evolutionary algorithms.— 1 ^ Therefore a cross-check of the proposed groundstate structures 
with alternative methods is always necessary to ascertain the results. With systematic and 
accurate first-principles calculations using evolutionary searching method combined with 
particle swarm optimization algorithm, 29 we confirm in this work that Cs-IV atomic phase 
with a space group of IAi/amd becomes the most stable phase beyond 490 GPa, but a novel 
structure Fddd is also found to be degenerate with it over a wide range of pressure. These 
two metallic phases persist up to 2.3 TPa, where exotic multi-atomic structures with low 
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symmetry becomes competitive. Our calculation also indicates that dense hydrogen has 
many very broad and flat basins on the energy surface, which has never been noticed before. 
This flatness of the energy variation not only leads to great anharmonicity in zero-point 
vibrations of lattice, but also blurs the boundary of some crystal structures. Furthermore, 
the transition paths between high pressure phases are modeled with nudged elastic band 
(NEB) method, and the associated enthalpy barriers and energy variations are analyzed. In 
the next section we will present the methods for total energy and NEB calculations. The 
results and discussion are given in Sec. IHIl followed by a summary in Sec. [TV] 

II. METHODOLOGY 

A. Total energy calculation 

The total energy calculations were performed with DFT as implemented in the Vi- 
enna Ab-initio Simulation Package (VASP).— The electronic structure was described with 
all-electron like projector augmented-wave (PAW) potential,— 1 ^ and the Perdew-Burke- 
Ernzerhof (PBE) exchange-correlation functional was used.— A hard version of the PAW 
potential that is specially designed for high pressure applications was employed. The k-point 
sets were generated separately for each unit cell encountered during the procedure, and a 
high quality Brillouin zone sampling with a grid ofl5xl5xl5 were found to be suffi- 
cient for structure optimization. When re-calculate the enthalpy curves, we used a denser 
k-point mesh that can generate at least 1500 irreducible points. The residual Pulay stress 
was removed by increasing the kinetic energy cutoff of the plane wave basis set to 900 eV, 
which was confirmed by observing the vanishment of the difference between the Hellmann- 
Feynman pressure and that computed from the energy curve with P = —dE/dV. Increase 
the energy cutoff to 1200 eV led to a pressure change less than 0.1 GPa. A justification of 
the methods, especially DFT and the employed pseudopotential, is given in appendix. With 
this parameter setting, the structural features are estimated to be converged to better than 
0.2%, the convergence of total energy is to within 3 meV per proton, and the relative energy 
difference between structures is to within 1 meV per proton. 

By minimizing the enthalpy, we carried out extensive and systematic structure searches 
of molecular and atomic phases at pressures up to 3.5 TPa. All structures were fully relaxed 
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at fixed volumes with a force tolerance of 0.1 meVA -1 . The pressure was directly computed 
with the Hellmann-Feynman theorem, which then was used to calculate the enthalpy. In 
addition to diatomic molecular phases, many high-pressure candidates of atomic phase were 
considered and their enthalpies were calculated, including cubic structures (SC, BCC, and 
FCC) and their low-symmetry distortions (/3-Po, /3-Hg, In-I, and In-II), hexagonal structures 
(HCP and w-Ti), diamond structure, IAi/amd (/3-Sn and Cs-IV) phases, and so on.— Details 
are presented in the following subsections. 



B. Structure search 

The search for high pressure structures of dense hydrogen was mainly carried out with 
the particle swarm optimization (PSO) technique 29 within an evolutionary scheme that 
combined with first-principles total energy calculations using VASP. PSO is an efficient 
approach of evolutionary methodology but quite different from genetic algorithm (GA). 
In particular the major evolution operations of crossover and mutation in GA have been 
avoided. PSO has been verified to perform well on many optimization problems.— >2L2£ In 
this work, stable and meta-stable structures containing up to 24 atoms per unit cell at 
pressures up to 3.5 TPa were automatically explored and generated by CALYPSO code,— 
which implements PSO algorithm. With a population size of 30 in each generation and a 
total allowed number of generations of 30, it is sufficient to ensure the convergence of the 
structure search. 

As a complement to the CALYPSO search, in order to better understand the relation- 
ship between symmetry and the structure stability, as well as to track the evolution of 
structures with pressure, we also performed local structure searches manually. It was based 
on the already known information about stable and meta-stable phases at relatively low 
pressures.— 1 ^ Some molecular and atomic phases which were not yet fully investigated pre- 
viously were added into the structure library of search. The atomic phases were selected from 
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by consulting previous theoretical predictions to pick out the most likely candidates, 
while putting high weight on low-symmetric structures. 

The initial structures for optimization were generated from the search library by drifting 
the ionic positions and distorting the lattice shape randomly. Interpolated configurations 
between the most stable phases were also taken into account. Although most of the phases 



in the search library are locally stable, with large enough distortions new structures can 
always be produced. To eliminate a possible sensitive dependence of the final results on 
the initial configurations, the ionic relaxation procedure employed algorithms of both quasi- 
Newton and conjugate-gradient method alternatively. The derived structures were then 
fully relaxed to the energy minima at constant volumes without any symmetry restrictions. 
After the ionic optimization process converged, we calculated the pressure and then the 
corresponding enthalpy. At each pressure, the search process was repeated several times 
(each iteration involves all structures in the search library that was updated accordingly) 
to confirm the results, and in total 4 typical pressure values (around 0.5, 1, 1.5 and 2TPa, 
respectively) were explored. The whole enthalpy curves of the low-lying phases were then 
computed. 

C. NEB calculation 

Structural phase transition paths and the associated energy barriers were modeled using 
the nudged elastic band (NEB) method^ as implemented in the VASP code, which searches 
for the minimum energy path (MEP) by moving a chain of ionic configurations or images 
that bridges the initial and the target structures. Tangential springs were introduced to 
keep the images being equidistant during the relaxation. The potential energy maximum 
along the MEP is the saddle point energy which gives the activation energy barrier. 

In calculations, the supercell method was employed, which in most cases contains 12 
atoms, and in some special situations a cell with 24 atoms was also used. To model the 
transition path with a variable cell volume and shape using NEB technique, the continuous 
MEP was constructed by linking the images with springs only in a pre-aligned supercell 
(i.e., the images due to periodic boundary conditions were discarded). When computing the 
transition path from high-pressure dense phases to a diatomic molecular phase at GPa, 
the molecular structure was generated by distorting and then relaxing the corresponding 
dense hydrogen phase respectively. This strategy simplifies the NEB calculations greatly 
since in this way the initial path generated through a linear interpolation of the initial and 
the target configurations gives a good approximation to the final transition path. Then a 
minimization of the whole system was carried out to trace out the MEP. In all calculations 
at least five images (seven if includes the two end images) were employed, which in most 
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cases is sufficient to give an acceptable resolution to the discrete representation of the MEP. 
D. Calculation of phonons in harmonic approximation 

Zero-point (ZP) vibrations of protons were neglected during the structure search and 
optimization procedure, but with a subsequent estimation of its impact on the relative sta- 
bilities using harmonic ZP energy (ZPE) calculations. The full phonon spectra in a harmonic 
approximation were calculated with the small-displacement method as implemented in the 
PHON code.— Big enough supercells containing more than 64 atoms were used. In the as- 
sociated DFT calculations, a Brillouin zone sampling mesh of 20 x 20 x 20, a kinetic energy 
cutoff of 1000 eV, and a very dense support augmentation charge grid that is required for an 
accurate force calculation were used. This setup gives a convergence in the ZPE better than 
2 meV per proton. The magnitude of the small displacement was slightly varied to check the 
numerical stability of the calculated force constant matrices. The ZPE was estimated from 
the phonon density of states g(oo) by J g(u)hu/2 du, and the ZP mean square displacement 
(MSD) was calculated by J g(u)h/ (2Mu) du, where M is the mass of hydrogen atom. 

III. RESULTS AND DISCUSSION 
A. Static structure and enthalpy 

1. Ground-state of dense hydrogen 

The calculated enthalpy difference with respect to a reference state (virtually defined by 
its enthalpy as H = — 3.985 + 0.216P 55 eV/proton, P in a unit of GPa) 37 of the most stable 
phases are shown in FigJTJ Overall, our results agree very well with previous theoretical 
predications: diatomic C2/c is the groundstate at 110 GPa and transforms into Cmca-12 
at about 255 GPa;— ^ at 370 GPa, Cmca overtakes slightly,— but it soon degenerates into 
Cmca-12, which is stable until an atomic phase- Cs-IV with space group I4i/ 'amd^ and its 
distortion Fddd- becomes the groundstate at 495 GPa; beyond 2.3 TPa, the most favored 
phase is shifted to trigonal i?3m in a static lattice approximation.— 

The main transition points are indicated by arrows in FigJIJ Point B is an intersec- 
tion of the molecular phases and the low-lying atomic Cs-IV and Fddd phases. Namely, 
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FIG. 1. (color online) Enthalpy difference per proton as a function of pressure with respect 
to a reference state. Arrows indicate the main phase transitions of the ground state and the 
shadow marks the region where the meta-stable low-symmetry multi-atomic phases lie in. Typi- 
cal curves are labeled as follows: open square-C2/c, open circle-Cmco-12, open triangle-Cmca, 
filled st&r-Pmmn, filled square-C2/m(2), filled up-triangle— C2/m(l), filled circle-/3-Hg, half-filled 
square-Cs-IV, half-filled triangle-Fdriri, filled down-triangle-R3m, crossed rhombus-i?3m, crossed 
pentagram-P63/mmc, crossed triangle-BCC, crossed circle-FCC, and crossed square-HCP. 

it is a pressure- induced dissociation point. This is in agreement with previous DFT 
calculations,— 1 ^ 1 ^ except for a newly discovered degenerate phase Fddd that was not 
detected before. Fddd is an orthorhombic variant of the diamond phase, but is also close to 
being a distortion of Cs-IV locally. Beyond the transition point C, there are many structures 
with closely competitive enthalpy, reflecting the frustration among competing factors. In 
this regime, a small change in interactions will tip the balance from favoring one phase into 
another. Although R3m seems to have the lowest static lattice enthalpy, harmonic phonon 
calculations indicated that the zero point vibrations of lattice might make R3m be more 
favored.— It is interesting to note that both Cs-IV and Fddd have 4 nearest neighbor (NN) 
atoms, and R3m has 6 as its first coordination number (CN). In contrast, R3m apparently 
has just 2 NNs, but since the distance from its first NN shell to the next one is very short 
(0.06 A at 3 TPa), so that it in fact has 6 atoms in its first coordination shell. From this point 
of view, the evolution of the groundstate structures of dense hydrogen under compression is 
evident: with increasing density, the structure evolves from with 1 NN (diatomic molecular 
phases) to that with 4 NNs (atomic Cs-IV and Fddd), and then to that with 6 NNs (R3m 
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FIG. 2. (color online) Structure of (a) C2/m(2), (b) C2/m(l), and (c) /3-Hg. A charge density iso- 
surface also illustrates in (c). The two C2/m phases (both in Immm representation and projected 
onto the same plane) are different only in the connection manner of the covalent bonds, and /3-Hg 
is derived from BCC structure with covalent bonds formed along the shortest c direction. 

or R3m), finally transforms to cubic structures with 8 (BCC) or 12 (FCC) NNs. 19 



2. Metastable multi-atomic molecular crystals 

Other monatomic phases with high symmetry, such as BCC, FCC, and HCP, and their 
distortions, do not approach the ground-state line before 3.5 TPa (the simple cubic (SC) 
phase is always 50 meV higher than others and not shown here). One interesting phase is 
/3-Hg, which caught no attention before. Its enthalpy is already low enough at 600 GPa 
as shown in FigfTJ and approaches gradually to the groundstate line all the way beyond 
3 TPa. It might coexist with other meta-stable phases such as Pmmn and C2/m over a 
wide range of pressure. In conventional sense it should be a monatomic phase distorted 
from BCC structure with a shorter lattice length in the c direction.— But for hydrogen at 
pressures beyond 600 GPa this lattice length is already short enough so that covalent bonds 



are established along this direction and it becomes a chained molecular state (see Fig 2(c)). 
In fact it is exactly the formation of this kind of multi-atomic molecular bonds that stabilizes 
the structure, due to a combined effect of pressure induced excitation of the Is electrons 
(which weakens the diatomic covalent bond) and the tendency to overlap electronic orbitals 
among neighboring atoms. It can be seen much clearly by comparing the enthalpy with 
that of Cmca (diatomic bonds), Pmmn (triatomic bonds, see below), and C2/m (molecular 
chains). Specifically, /3-Hg constitutes linear chains of H 2 bonds, whereas C2/m consists of 
chains of H3 clusters. C2/m(l) and C2/m(2) belong to the same space group. The only 
difference is in the linkage pattern of the bonds and the chain orientation. In addition, 
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FIG. 3. (color online) Charge density on [010] plane (a) and the atomic structure (b) of a 2x1x2 
supercell of Pmmn phase at 1104 GPa, in which the H3 clusters are evident. In (b): atoms on 
alternating layers are distinguished by different color/grey scale, and the dashed lines indicate the 
connection manner of covalent bonds in a chained molecular phase C2/m(2). 



C2/m(l) has a shorter interchain distance and trends to form a bond network, whereas 
C2/m{2) has more localized bonds and thus its structure is much distinct and well defined. 
Figures EH3] and table U in appendix provide detailed structural information about these 
phases. It is worthwhile to point out that the C2/m phase reported in RefJl9l corresponds 
to the C2/ra(l) here, and C2/ra(2) has a lower enthalpy at high pressures. The existence 
of this isomorph (and others alike) indicates the complex nature of the structure of dense 
hydrogen. 

Although Pmmn is a triatomic phase (see FigJH]) and C2/m(2) is a chained molecular 



structure (see Fig, 2 (a) ), they are closely related. Both phases can be drawn on an orthorhom- 
bic lattice, and their structural relationship manifests clearly in an Immm representation. 



By comparing Figs. 3(b) and 2(a) (and also their charge densities), it is easy to find the 
similarity of these two structures: by compressing the former phase along its c direction to 



create new bonds along the dashed lines as indicated in figure 3(b) meanwhile shifting the 
relative position of the alternating layers slightly along this direction, one gets the latter 
phase. This is actually how the structural transition from Pmmn to C2/m(2) takes place 
at 1.9 TPa. It is reasonable since compression reduces inter-distance among H3 clusters, and 
thus could establish chains of covalent bonds along certain direction if the strain tensor is 
anisotropic. 

It is worthwhile to note that at the intermediate pressure range the NN number of C2/m 
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phases is about 3, and both Pmmn and /3-Hg have a NN number of 2. They all are multi- 
atomic molecular phases. Furthermore, distorting some simple high-symmetric structures 
can make them continuously relax to configurations with a low CN, and the resultant en- 
thalpy lies in the shaded range as shown in Fig{TJ The above mentioned multi-atomic 
molecular phases are members of these low-symmetric structures. This connection provides 
a plausible answer to the questions of why structures with a CN of 2 or 3 absent from 
the groundstates, and why dense hydrogen dissociates into orthorhombic atomic Cs-IV in- 
stead of simple high-symmetric monatomic phases such as BCC, FCC, or HCP. As it is well 
known, hydrogen has just one electron outside the nucleus. This makes the electron cloud be 
compact and tightly attached to the nuclei even at high pressures. On the other hand, at a 
given density, high-symmetric monatomic phases always have a greater interatomic distance, 
making it difficult for hydrogen to share electron with its neighbors. In this sense, distortion 
of the structure to reduce interatomic distance between some atoms can effectively facilitate 
wavefunction overlapping and lower the energy. That is the reason why dense hydrogen 
does not dissociate into simple monatomic phases directly. From another point of view, 
diatomic covalent bond has been weakened greatly by compression at high pressures, and 
thus prevents a further transition of these low-symmetric phases into a diatomic molecular 
configuration via mechanism- for example, the Peierls instability^- from happening when 
the pressure is higher than 500 GPa. However, because each hydrogen has only one electron, 
it is not easy to form stable multi-atomic (triatomic or molecular chain) bonds. Even with 
the aid of a small portion of p electrons that are excited from s orbitals by compression 
(this increases the anisotropy of the electron cloud and thus favors multi-atomic bonds), the 
multi-atomic molecular phases are just be meta-stable, and the balance can be easily tipped 
down towards the atomic phases with a low CN, which takes advantages of all competing 
factors. That explains why no groundstate of dense hydrogen has a NN number of 2 or 3. 



3. Electronic structure 



Figure 3(a) shows the calculated charge density of Pmmn at about 1 TPa. The tri- 
atomic molecules are evident. Note that the environmental charge density is nonzero and 
the multi-atomic molecules are in fact submerging in a sea of electrons that mediates metallic 
interactions. Figure H] plots the electron localization function^^ of C2/m(2) at 2.7 TPa, 
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FIG. 4. (color online) Electron localization function of the C2/m(2) phase (in a tetragonal Immm 
representation) at about 2.7 TPa. Formation of the molecular chains is evident. 

where the weak covalent bonds along the molecular chains built up of H3 clusters are distinct. 
It is important to point out that although C2/m and /3-Hg are chained molecular phases, they 
are natural conductors. We do not need to perform a band structure calculation to confirm 
this, since these phases have an odd number of electrons in their primitive cell, therefore 
there always have a half-filled band. Although Pmmn has 6 electrons in its primitive cell, 
a DFT calculation with GGA showed that it is metallic when pressure is beyond 600 GPa. 
This is different from lithium or sodium where intermediate compression induces localization 
of electrons at interstitial regions and leads to metal-semiconductor-metal transition.— 

The electronic density of states (DOS) of C2/c, Cmca-12, Cs-IV, and C2/m{2) phases at 
selected pressures are shown in Figj5j respectively. Being consistent with previous studies, 
the molecular hydrogen at 299 GPa does not show metallic characteristics, and there is an 
energy gap presented at the Fermi level.— At a pressure of 490 GPa in a diatomic molecular 
Cmca-12 phase, the gap already closes up and the material becomes metallic. Here we 
didnot attempt to determine the precise insulator-metal transition pressure in the molecular 
phase but instead focused on the general trend of the variation of DOS with pressure. 
It is interesting to note that the DOS dips down at the Fermi level in Cmca-12, clearly 
showing that the closure of the gap is due to band overlapping. This feature disappears 
after transition into the atomic Cs-IV at the same pressure, it also absents from C2/m(2) at 
about 2.7 TPa in spite of covalent bonds presenting in this meta-stable phase. Overall, with 
increase of the pressure, the localized covalent states in diatomic molecular phases become 
dispersive and extend towards both the high and low energy ends, which then closes up the 
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FIG. 5. Electronic density of states of dense hydrogen at high pressures. The dash-dotted vertical 
line indicates the Fermi level. Notice that C2/m(2) at a pressure of 2.7 TPa has typical features 
of simple metals. 

gap by band overlapping at the Fermi level. On the other hand, Cs-IV and C2/m(2) show 
typical characteristics of simple metals in their DOS. Especially, the DOS of C2/ra(2) has 
already been dispersed greatly by compression and becomes flat and featureless over a wide 
range of energy. Other phases of dense hydrogen at ultra-high pressures are similar and it 
is unnecessary to discuss them separately. 



B. Effects of zero-point motion of protons 

1. Harmonic ZPE and its failure 

The contribution of zero-point vibrations of protons to energy and enthalpy can be taken 
into account within quasi-harmonic approximation. Since in a static lattice approximation 
our calculated groundstates and the relative stable order of high pressure phases of dense 
hydrogen are almost the same as that reported in Ref , hoi, and they had performed a detailed 
analysis of the zero point energy (ZPE) in harmonic approximation, thus it is not necessary 
to repeat the discussion. Here we only illustrate the magnitude and the possible consequence 
of harmonic ZPE using meta-stable multi-atomic molecular phases as examples. As Figj6] 
shows, inclusion of harmonic ZPE slightly changes the relative stability of these phases. The 
magnitude of ZPE is about 0.32 eV per proton at 1 TPa, and increases to 0.42 eV per proton 
at 2 TPa. For highly compressed atomic phases of hydrogen, this treatment is inappropriate 
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FIG. 6. (color online) Change of the relative stability of meta-stable dense hydrogen due to 
contribution from harmonic ZPE, where the solid (dased) lines indicate with (without) ZPE, AH 
stands for the enthalpy difference with respect to a predesignated reference state, and the filled 
circles, squares, triangles, and half-filled squares denote C2/m(2), C2/m(l), /3-Hg, and Pmmn 
phases, respectively. 

because of the strong anharmonic effects.— _ — However, for multi-atomic molecular phases, 
the anharmonic effects are also remarkable and the shape of the potential well around the 
equilibrium position of each proton is far from being a quadratic form (see below), which 
undermines the justification for harmonic approximation. In this sense the relative stability 
of dense hydrogen cannot be faithfully determined by harmonic ZPE, because the enthalpy 
difference between these phases is too small, and at the same time it is almost impossible 
for the harmonic approximation to have a precision of within several percents in the case of 
dense hydrogen. 

On the other hand, most of the energy barriers that separate different phases of dense 
hydrogen are much smaller than the magnitude of harmonic ZPE. In this situation, the 
classical notation of structural phase might be ill-defined, because it is easy for the system 
to overcome the barriers and travel freely from one phase into another driven by ZP mo- 
tions. Such kind of quantum fluctuation between structures invalidates not only harmonic 
approximation of lattice dynamics, but also some restricted quantum Monte Carlo (QMC) 
treatments.— 1 ^ To tackle this problem quantitatively, a full quantum treatment of protons 
on the same footing as electrons is required,— 1 ^ which however, is difficult within the DFT 
framework based on Born-Oppenheimer approximation. 
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2. Anharmonic zero-point motion 

Although in DFT it is difficult to carry out a quantitative analysis of the anharmonic 
ZP effects if nuclei are treated as classical particles, an insightful perception of the quantum 
motion of protons still can be obtained by inspecting the energy surface or landscape closely. 
Figure [7] shows a series of section of the energy surface that cut along a transition path from 
FCC to /3-Hg at different densities. That is to say, changing the lattice length c of the 
/3-Hg structure while adjusting the a and b accordingly at a fixed density with r s equaling 
1.12 (~lTPa), 1.04 (~1.5TPa), and 1.01 (~2TPa), respectively. Here the dimensionless 
parameter r s is defined as the radius of a sphere which encloses on the average one electron 
in a unit of the Bohr radius.— - — A value of r s — 1 corresponds to a physical compression 
ratio of about 30 for hydrogen. The energy surface obtained in this way is approximation 
free, except those that already introduced in the standard DFT formalism. Particles move 
on these surfaces at K as zero point motions. 

It is easy to find from these surfaces that the harmonic approximation breaks down com- 
pletely for monatomic phases. There is even no potential well can be defined between BCC 
and FCC (and between BCC and /3-Hg as well). The flatness of the energy surface unveils 
an important origination of the unusually large anharmonic effects in monatomic phases 
of hydrogen observed in QMC calculations (in addition to the light mass of protons):— - — 
having such a flat energy surface, the crystal might melt even at zero Kelvin. 

On the other hand, the harmonic ZP root mean square displacement (rMSD) gives a 
typical value of 0.2 A for /3-Hg phase.— It characterizes a typical ZP motion size for protons 
in a chained molecular phase. Making use of this value and the potential well width given in 
FigJTJ we estimated that at a pressure of 1 TPa the ZP motions of protons are confined within 
the potential well. But with increasing of the pressure, ZP vibrational size becomes com- 
parable with the potential well width at 1.5 TPa, and exceeds the latter at 2 TPa. Namely, 
beyond that pressure /3-Hg might merge into FCC and BCC phases driven by ZP motions. 
This picture can be understood by analogy with an effective particle that moves on a similar 
energy surface (a half-infinite potential well), as shown in the inset of FigfTl In that case a 
bound state (i.e., a well-defined crystal) exists only when 8MVx 2 > ir 2 h 2 where x(V) is the 
potential width(depth) and M is the particle mass.— When no bound state is available, the 
particle moves around freely, but a phase cancellation between the incident and reflected 
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FIG. 7. (color online) Energy variation with the change of the c axial length of the /3-Hg structure, 
which passes through FCC, BCC, and /3-Hg successively. Inset: a simplified half-infinite potential 
well that models the energy surface, as well as the corresponding probability density profile of an 
effective particle moving in it. 

waves make it have a relative high probability within the potential well, as the probability 
density profile in the inset of FigJT] illustrates. 

Alternatively, this phenomenon can also be understood intuitively via the path integral 
formalism of quantum statistics theory. The partition function of a system consisting of 
distinguishable particles (here hydrogen atoms) can be written as 51 

(1) 



Z = J p(R, R; P) dR, 
with the diagonal density matrix given by 



p(R,R;(3) = p°(R,R;(3)(ex P 



V{Rt) dt 



(2) 



BW 



Here p° is the free particle density matrix, /3 denotes (&bT) _1 where T is the temperature, 
and the potential contribution (the () bw term) is given by averaging over Brownian motions 
along closed paths that weighted by potential energy. The total energy with ZP contribution 
included is thus 

Eo = - km —-. 3 

^->oo dp 

In cases where structural phases are separated by high enough barriers, the energy E is 
locally defined. Namely, all Brownian motions in Eq.(J2]) are effectively confined to a limited 
phase space which characterizes the structure, and therefore one can compare the value of 
Eq to determine the relative stability between different phases. 
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However, when the energy barrier is low or even no barrier at all, Eq of one phase contains 
contribution of the Brownian paths winding across other phases, and the definition of the 
energy of that phase becomes meaningless. In this situation one could instead analyze the 
probability for a Brownian motion to fall into a specific region in the phase space. Taking 
/3-Hg for example, the probability for this phase to appear is given by 

_ In m gj g) dR 
V J n p(R,R;P)dR' U 

where fl m is the domain of the phase space defined by the structure of /3-Hg, i.e., an analogue 
of the width of the potential trap as shown in FigJTJ and Q is the whole phase space domain 
that is accessible to all Brownian motions. Because the potential well is attractive, the paths 
that belong to/or pass through the potential well (fl m ) always have a higher weight than 
others. Therefore rj has a non-zero value at low enough temperatures. With the temperature 
decreases further, Q gradually shrinks backwards to Q m . At zero Kelvin, if Q reduces to 
Q m exactly, then the phase of /3-Hg is a well defined classical structure. Otherwise there are 
quantum fluctuations and /3-Hg exists only instantaneously with a probability of rj. Note this 
kind of quantum behavior that blurs the structural boundary of phases had been noticed 
in path integral simulations,— 1 ^ which treated the quantum nature of protons. Here we 
considered only meta-stable phases, similar argument (though not identical) can be applied 
to the groundstate Cs-IV and its distortions. 

C. Energy barrier of phase transition 

From above discussion, we knew that the shape of the energy surface and the barrier 
that separates different phases are crucial content to understand the high pressure behavior 
of dense hydrogen comprehensively. The transition path shown in FigJT] is a special case, 
which reveals that no energy barrier presents between chained /3-Hg and cubic FCC or BCC 
phases. Furthermore, the potential well around /3-Hg becomes shallower and shallower with 
increasing pressure, indicating the enhancement of the stability of high-symmetric phases 
and the weakening of the tendency to form low coordinated structures. It should be noted 
that this path is along a preassigned route. We can do this because the transition path is 
fixed completely by symmetry and structure analysis. However, this is not always the case, 
and a general technique such as NEB has to be employed in order to map out the transition 
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FIG. 8. (color online) Energy and enthalpy variations along the phase transition path from Pmmn 
to Cs-IV at 1 TPa. Note that Cs-IV has a freedom to drift far away from its ideal position. 

path. 

Figure M illustrates the transition path from meta-stable Pmmn to the groundstate Cs-IV 
phase at lTPa calculated with NEB. Both variations of energy and enthalpy are shown. A 
distinct barrier (0.035 eV/proton in enthalpy and 0.055 eV/proton in energy from the Cs-IV 
side) was clearly obtained. From the energy variation, we can see that both of the end phases 
have distortions with a small energy change. In particular, the groundstate Cs-IV has various 
variants which possess almost the same energy or enthalpy. The average drift distance of 
atoms between these variants can be as large as Ad = 0.39 A, compared with the NN distance 
of 0.92 A in this structure. This implies that it might be quite common for broad and flat 
basins to present at some locally stable phases (even the groundstate) of dense hydrogen, 
which was never noticed before. Furthermore, the flatness of the basins implies that the ZPE 
should be much smaller than that predicted by harmonic approximation. In addition, since 
an enthalpy barrier of 0.035 eV/proton from Cs-IV to the first meta-stable phase Pmmn 
corresponds to a temperature scale of 400 K, it is highly possible that dense hydrogen is in 
a solid state at this pressure range at room temperature. Of course, considering that Cs-IV 
itself has many variants with negligible energy change, it is also possible that it melts locally 
driven by ZP motions. If this is the case, then a liquid-liquid transition could be expected 
when temperature is increased. 

Being analogous to Figd a well-defined barrier is absent between the meta-stable Pmmn 
and C2/m(2) at 2.8 TPa, as shown in Figj9j There are three noticeable characteristics 
presented in this path: (i) a more stable R3m phase manifests itself in the transition path; 
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FIG. 9. (color online) Energy and enthalpy variations along the phase transition path from Pmmn 
to C2/m(2) at 2.8 TPa. Note that the energy variation does not comply with that of enthalpy, and 
a more stable phase of R3m is unveiled. 

(ii) a flat energy (or enthalpy) surface appears again; and (iii) the energy variation is inverse 
to that of enthalpy. The first point shows that it is possible to find out more stable structures 
by investigating the transition path between high-lying meta-stable phases. This is also 
helpful for understanding the physical mechanism of structure stability. The second point 
reflects the frustration of different competing factors, and implies that it is ineffective to 
optimize the structure by conventional relaxation algorithms because the forces become 
too small to evolve the geometry within these flat regions. The final point clearly implies 
that C2/m(2) is stabilized by the term of PV, and should have great imaginary phonon 
frequencies in harmonic approximation because the lattice dynamical matrix is evaluated on 
the energy surface. However, such imaginary modes do not necessarily mean that the phase 
is locally unstable in thermodynamics if one takes the PV contribution into account. 

It is interesting to investigate the transition path from diatomic molecular Cmca- 12 to 
atomic Cs-IV at the dissociation pressure of 495 GPa. As Fig(T0] shows, a high enthalpy 
barrier of AH = 0.038 eV/proton presents. This will introduce hysteresis and lift the ap- 
parent dissociation pressure at low temperatures. Since the volume difference per atom 
between Cmca-12 and Cs-IV is about AV^ = 0.03 A 3 , the possible pressure increase is thus 
AP = AH/ AV~ 200 GPa to the first order of correction.— That is to say, the DFT pre- 
dicted dissociation pressure of diatomic hydrogen should be at about 700 GPa if taking the 
kinetic effect of enthalpy barrier into account. Other factors that might also have some 
impacts include: (i) ZP motion of protons, which favors atomic phases and therefore should 
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FIG. 10. Enthalpy variation of dense hydrogen along the transition path from diatomic molecular 
Cmca-12 to atomic Cs-IV phase at 495 GPa, where a low-symmetric filamentary phase and Fddd 
phase (a distorted diamond structure) were observed. 



decrease the dissociation pressure,— dMl£2. but its precise value is hard to evaluate at present; 
(ii) intrinsic DFT error which favors homogeneous distribution of electronic density, thus 
also underestimates the dissociation pressure about 50 GPa;^ (iii) other transition path with 
lower enthalpy barrier, we cannot exclude this possibility completely because NEB is in fact 
a local optimization algorithm, and it explores a limited phase space and thus depends on 
the initial guess of the transition path in some degree. Nevertheless, after taking all of these 
factors into account the dissociation pressure should be less than 750 GPa, or 550 GPa if no 
any hysteresis effect presents. 

Another interesting phenomenon unveiled in Fig{10] is that there is no energy barrier 
between the degenerate groundstates of Cs-IV and Fddd. In other words, dense hydrogen 
fluctuates between these two structures driven by ZP motions. Note that although Fddd 
can be viewed as a distortion of Cs-IV structure, from a crystallographic point of view, it 
is in fact an orthorhombic variant of the diamond phase. Namely, the quantum structural 
fluctuation in the groundstate of dense hydrogen is much more drastic than what FigJH] has 
been implied. Also an anisotropic low-symmetry filamentary phase manifests itself on the 
transition path, which is very similar to the previously proposed low-symmetric structures 
of dense hydrogen.— ^£~— It is separated from other phases by an energy barrier and might 
be meta-stable. 
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FIG. 11. (color online) Energy variation of a series of dense hydrogen phases along the respective 
transition path to a diatomic molecular structure at zero pressure. No energy barrier can be 
detected. 

D. Stability of dense hydrogen at zero pressure 

The stability of atomic phases of dense hydrogen at zero pressure is intriguing because 
it is a possible room temperature superconductor.— ~— Although early primary perturba- 
tion calculations on the structural energy of hydrogen within a static-lattice approximation 
with the effective electron-ion interaction expanded to fourth order suggested that a highly 
anisotropic structure might be stable,— 1 ^ later more careful analyses dismissed this pro- 
posal and found that isotropic phases are more favored.— iMiOl The energetic stability of 
these recently proposed new groundstates and low-lying meta-stable phases^^ 1 ^ at GPa, 
however, has not been studied yet. Namely, whether these high pressure phases of hydrogen 
is quenchable or not is still unknown. 

To answer this question, we investigated the transition path and the associated energy 
barrier of these dense phases to diatomic molecular structures at GPa with NEB calcu- 
lations. The initial dense structures were prepared by carefully relaxing the configurations 
from high pressures with constraints applied. The corresponding diatomic molecular phases 
were then produced by distorting and relaxing the respective dense structures. Typical re- 
sults are illustrated in FigJPTl The conclusion is that all of the low-lying states of dense 
hydrogen known so far, including isotropic atomic phases and low coordinated tri-atomic or 
chained molecular phases, have no transition barriers can be detected along the transforma- 
tion path. In particular, the exotic filamentary structure observed in FigJlO] has no energy 
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FIG. 12. Calculated equation of state at OK of dense hydrogen up to 3.5 TPa compared with the 
extrapolation of the principal Hugoniot. Inset shows the comparison with the Vinet EOS (dashed 
line) that fitted to experimental data up to 119 GPa at low pressures. 

barrier, too. That is to say, there is nothing that can prevent these phases from sponta- 
neously decaying to diatomic molecular structures. This fact that dense hydrogen might 
be unquenchable is a direct consequence of the strong tendency of hydrogen to pair at low 
pressures (Pmmn and C2/m are also unstable at OGPa, which didnot include in FigfTT]). 
This result is obtained with DFT in PBE approximation. This exchange-correlation energy 
functional constructed on a local or semi-local approximation to the homogeneous electron 
gas value overestimates the stability of a metallic phase slightly, but cannot eliminate the 
whole energy barrier completely if there were one. Also, an insightful analysis of the intrin- 
sic error in DFT— ^ implies that it seems unlikely that an additional energy barrier can be 
predicted by an exact theory of DFT when no any barrier can be detected with the current 
version of DFT. That is to say, even in a level of the exact many-body quantum theory, 
dense hydrogen might still be unquenchable for these already known structures. 

E. Equation of state at zero Kelvin 

Equation of state (EOS: the pressure-volume relation here) is an important content to 
understand the response of a material to compressions, which reflects not only the phase 
transitions driven by pressure, but also the evolution of interactions with the change of 
coordination environment. In experiment, x-ray diffraction measurements of the structure 
of single- crystal molecular hydrogen had been performed at pressures up to 109 GPa for H 2 
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and 119 GPa for D2.— From these measurements a high-pressure EOS was deduced,— which 
is significantly more compressible than that early constrained from lower-pressure data to 
42 GPa.- This deviation was interpreted as gradual effect of orientation of the H 2 molecular 
axis within phase I.— At higher pressures, hydrogen transforms into atomic Cs-IV phase, 
where both interatomic interactions and crystal structure are changed. Therefore there is 
enough reason to suspect that this experimental EOS cannot be extrapolated to higher 
pressure range. 

With the information of the calculated groundstate structures, the whole pressure- volume 
curve of dense hydrogen at zero Kelvin was calculated up to 3.5 TPa with DFT. Figure H~2l 
shows this curve by comparing with the extrapolation of the principal Hugoniot. It is 
worthwhile to note that although there is a little difference between R3m and R3m in the 
enthalpy and it is hard to determine which one is the true groundstate at ultra-high pressures, 
both phases have an almost identical P-V relation. Similarly, Fddd has the same variation 
of density with pressure as that of Cs-IV and its distortions, and thus didnot include in 
FigJT2l Along the groundstate line, the volume collapses due to phase transitions are small, 
it is about 2.59% when dissociates to Cs-IV phase and 1.22% between PQ^/m and C2/c. In 
other transitions it is just about 0.5%. 

The experimental data can be fitted to a Vinet function,— which is overall in good agree- 
ment with our DFT calculated data below 300 GPa, as the inset of Fig|T2] shows. However, 
there is a subtle discrepancy: the DFT overestimates the compressibility slightly when be- 
tween 50 and 150 GPa. This might be due to the fact that we didnot treat the quantum 
rotation of H2 molecules explicitly. This degree of freedom of motion should contribute 
energy and thus increase the internal pressure accordingly. At higher pressures, however, 
extrapolation of the fitted Vinet function overestimates the compressibility significantly. As 
Fig{13] illustrates, both of the Vinet function that fitted to original experimental data (Vinet 
300 K) and that fitted to data reduced to K without ZP contribution (Vinet K) deviate 
from the DFT results when beyond 500 GPa, and the difference can reach as high as 36% 
at a density of r s = 0.9. Note that the curves of "Vinet 300 K" and "Vinet OK" are almost 
identical within the studied pressure range, showing that both the ZP pressure and ther- 
mal pressure are relatively small in comparison. Furthermore, since in DFT calculations 
we didnot take ZP contribution into account and usually ZP motion contributes a positive 
pressure and thus should reduce the compressibility, the DFT EOS shown in Fig{13] is thus 
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FIG. 13. Comparison of the calculated pressure-volume curve with the Vinet EOS models that 
fitted to low-pressure data and a proposed EXP/P5 function within a pressure range up to 3.5 TPa. 

an upper estimate of the compressibility of dense hydrogen, and any model with a higher 
compressibility (e.g., the fitted Vinet EOS) is not allowed in physics. 

If ignored the small volume collapses at the phase transition points, the pressure- volume 
relation of the groundstates of dense hydrogen at OK calculated by DFT can be fitted 
excellently by a function of 

5 

P = Y[lO Anr ". (5) 

n=0 

We denote it as EXP/P5. When the pressure P is in a unit of GPa, the parameters are as 
follows: A = 1.0683, A x = 19.1824, A 2 = -36.3776, A 3 = 28.5165, A 4 = -10.6068, and 
A 5 = 1.5224. It should be noted that EXP/P5 function not only faithfully represents the 
overall variation of the EOS of dense hydrogen over a broad range of pressure up to 3.5 TPa, 
it also reproduces the low pressure data accurately, as demonstrated in the inset of FigJT2l 
With these properties, this function should have a good applicability for extrapolating the 
EOS of hydrogen to higher pressures beyond several TPa. 

IV. CONCLUSION 

High pressure phases of dense hydrogen have been extensively searched up to 24 atoms 
per unit cell and 3.5 TPa in pressure with density functional theory calculations. The results 
confirmed the previous conclusion that diatomic hydrogen dissociates into atomic Cs-IV at 
495 GPa and then transforms to R3m or R3m at about 2.3 TPa. Exotic high pressure 
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behaviors were also discovered that show dense hydrogen having numerous isomorphs with 
comparable enthalpy, and frustrated competitions among them lead to broad and flat basins 
on the energy surface. In particular, calculations showed that there is even no energy 
barrier to separate crystallographic different phases in some cases. Furthermore, the atomic 
groundstate Cs-IV has a freedom to distort and is degenerate with orthorhombic Fddd, 
which fluctuates forth and back to Cs-IV driven by ZP motions of protons. 

Within a wide range of pressure beyond 500 GPa, the first meta-stable structure is a tri- 
atomic Pmmn, therefore the stability of the closely related multi-atomic molecular phases 
was also analyzed. Electronic structure calculations indicated that they are metallic in na- 
ture, but with weak covalent bonds presented. The occurrence of these meta-stable phases 
is a direct consequence of the competition between the tendency to overlap the electron 
orbitals among neighboring atoms and the pressure-induced dissociation of molecular hy- 
drogen. The structural relationship between these exotic phases was investigated for a better 
understanding of the structural behavior of dense hydrogen. 

In addition to static-lattice calculations, ZP vibrations of lattice were also computed in the 
harmonic approximation and the magnitude of its contribution was estimated. Inspection 
of the energy surface (via transition path) dismissed the validity of this level approximation, 
and the anharmonicity of ZP motions was demonstrated by studying the variation of the 
energy surface with pressure, which not only reveals the weakening of the tendency to 
form covalent bonds under compression, but also illustrates the importance of quantum 
fluctuations in structure of dense hydrogen at high pressures. 

The general transition path and the associated energy barrier between different phases 
were calculated with NEB technique, which suggested that the dissociation pressure of 
diatomic molecular hydrogen might be deferred to at about 750 GPa due to a hysteresis effect 
of the energy barrier at low temperatures. The enthalpy barrier between the groundstate 
Cs-IV and the first meta-stable Pmmn phase at 1 TPa implies atomic hydrogen might be 
in a solid state at zero Kelvin, but it depends on the exact contribution of anharmonic ZP 
vibrations of protons. The calculated transition path also suggested that there are phases 
in dense hydrogen that are stabilized purely by PV term, which is unusual in the common 
sense. The meta-stability of dense phases of hydrogen at zero pressure were extensively 
studied. No energy barrier was detected between them and the diatomic molecular phase, 
implying dense hydrogen might be unquenchable. 
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The equation of state (i.e., the pressure- volume relation) of dense hydrogen up to 3.5 TPa 
was investigated. A new EOS function, namely EXP/P5, was proposed, which can repre- 
sent the calculated DFT data over the whole studied pressure range excellently. Since no 
thermal pressure and ZP contribution were included, this EOS is an upper estimation of the 
true compressibility, provided that the current understanding of the groundstate structure 
of dense hydrogen is correct. At low pressures, this EOS is in good agreement with the ex- 
perimental data measured by single crystal x-ray technique. Extrapolating the Vinet EOS 
model that fitted to these measured low-pressure data, however, drastically overestimates 
the compressibility when beyond 500 GPa, due to inappropriate counting of the interatomic 
repulsion under extreme compressions in this model. By the way, as isotopes of hydrogen, 
deuterium and tritium have almost the same electronic behavior. The main difference is 
that they have a heavier ionic mass, and thus the ZP effects should be less significant. But 
the overall picture of the physics is the same. 
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Appendix A: Justification to the density functional theory 

The sensitivity of the calculated results to the choice of exchange-correlation density 
functional had been checked in RefJl7l. We repeated the checking process and obtained 
similar conclusions. Furthermore, the validity of the density functional theory to ultra- 
high pressure physics can be understood easily. We know that the apparent failure of all 
local density approximation (LDA) based functionals relates directly to an abrupt variation 
of charge density. Hydrostatic compression reduces interatomic distance, which increases 
interactions between electrons and nearby nuclei. Covalent bonds are thus weakened and 
electronic wave-function spreads out due to environmental changes in atomic coordination 
usually. The direct consequence is thus a more smooth spatial distribution of the charge 
density, and a better performance of LDA-based functionals with an increase of the pressure 
(note that the peculiar phenomenon of pressure-induced localization of electrons observed 
in lithium2Z^i and sodium^ is absent from dense hydrogen). In addition, the s orbital 
electrons, the only one that is relevant in hydrogen, are always well described by LDA even 
at ambient conditions. 

On the other hand, although a reliable convergence of the DFT total energy to within 
3 meV per proton has been achieved, a concern that whether DFT can be applied to dense 
hydrogen beyond insulator-metal transition might still exist because it is well known that 
DFT underestimates band gap and then the transition pressure. It is a severe problem 
for electronic structure properties, but has limit impacts on the total energy and geometry 
features. Theoretically DFT only ensures the correct charge density and the total energy, 
rather than the quasi-particle levels that are introduced via Kohn-Sham ansatz. For example 
in an extreme case where the exact exchange-correlation functional were available, DFT 
would become rigorous and predict an exact charge density and total energy, whereas the 
band structure still cannot be generally guaranteed theoretically.— 1 ^ On the other hand, 
there are many examples in the applications of DFT to condensed matters where the band 
structure given by DFT is wrong but the energetics and atomic structure are still within an 
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FIG. 14. Variation of the tetragonal lattice parameters of /3-Hg under compression, where the 
shortest lattice vector c, which indicates the bond length of the molecular chain, is always longer 
than the ambient H2 value of 0.74 A. 

acceptable precision. Therefore the underestimation of the band gap of dense hydrogen by 
DFT does not present as a serious issue to what we concerned here. The intrinsic error in 
the total energy and structural features comes from the GGA approximation of the energy 
functional, and should be similar for insulator and metallic phases as long as the spatial 
variation of the charge density is similar. 

It also needs to point out that near the transition point from insulator to metallic state, 
where the band gap begins closing up and electrons jump abruptly from localized to de- 
localized orbitals, DFT might underestimate the transition pressure,- presumably due to 
self-interaction errors in exchange-correlation functional. But far away from this transi- 
tion region, DFT works well.- 1 ^ Since exact-exchange calculations predicted a metallization 
pressure of hydrogen at 400 GPa,— much lower than the pressure range interested here, we 
estimate that the influence on the total energy and structure features due to the errors in 
derealization of electrons occurred at low pressures should be small. 

Appendix B: Justification to the pseudopotential 

We used a PAW pseudopotential that is harder than the standard one, which is specially 
designed to account for the possible short interatomic distance under ultra-high compres- 
sions. The performance of this potential was checked by comparing the bond length and 
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binding energy of H2 dimer with those of other potentials and all-electron calculations. The 
low-lying structures and the corresponding transition pressures below 400 GPa as reported 



in Ref 



171 were also perfectly reproduced. This verified that the hard potential works cor- 



rectly. Its applicability to higher pressure range can be generally guaranteed as long as the 
shortest interatomic distance is greater than twofold of the outmost cutoff radius r m of the 
potential, which defines the atomic spheres (augmentation regions) where the pseudo wave 
function takes effects. This precondition ensures the correct wave-function shape near the 
atomic spheres. For comparison, r m =0.42 A for the hard potential and r m =0.58A for the 
standard potential of hydrogen. In our studied pressure range here, no low-lying structure 
has a shortest interatomic distance less than the ambient hydrogen molecule bond length 
of 0.74 A. For example, FigJT4l shows the variation of the lattice parameters of /3-Hg with 
pressure, where c indicates the shortest bond length and is always greater than 0.74 A. The 
overlapping between atomic spheres of dense hydrogen is less than that of H 2 molecule at 
ambient conditions. Therefore the application of the PAW potential of hydrogen to high 
pressures does not present as a difficult issue. It is worthwhile to point out that for am- 
bient H2, there is a little overlap of the pseudopotential regions on nearest neighboring 
hydrogen atoms. In principle these spheres shouldnot overlap, but a bit of overlap/softness- 
of-the-potential tradeoff is usually acceptable. Especially when spherical Bessel functions 
were used to construct the PAW potential, as implemented in VASP code, a large overlap 
between the atomic spheres is allowed.— 

An overall examination of the validity of the PAW potential is illustrated in FigfT5l 
where the VASP result is compared directly with that of the all-electron method calculated 
with WIEN2k code.™ Both calculations used the same PBE exchange-correlation energy 
functional. We can see from this figure that within a wide range of pressure up to 4 TPa (r s ~ 
0.9), both methods give almost the same energy difference between high-symmetric BCC 
and FCC structures (the maximum deviation is less than 0.4 me V per proton), reflecting 
the fact that the error introduced by PAW potential is insignificant. 



Appendix C: Structure of some meta-stable phases 

The structural information (including both lattice parameters and atomic coordinates) 
of some meta-stable multi-atomic molecular phases and the degenerate groundstate Fddd 
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FIG. 15. Comparison of the energy difference between FCC and BCC phases as a function of 
density calculated using PAW pseudopotential method with VASP from that of all-electron results 
of WIEN2k. 

at selected pressures are listed in table H 
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TABLE I. Structure of some low-lying phases. Only the fractional coordinates of symmetry in- 
equivalent atoms are reported. The numbers of atoms in a primitive cell are given. Immm is a 
coarse but robust (with greater tolerance) orthorhombic representation of C2/m structure. For 
/3-Hg, a tetragonal distortion of the BCC cubic cell is adopted. 



space group 


pressure 


lattice parameter 






atomic coordinates 


(j atoms 


(GPa) 


(A, °) 






(fractional) 


C2/m{\) 


650 


o=3.178 6=1.203 c= 


=2.395 


HI 


0.9999 0.0000 0.8300 


3 




a=7 =90.00 (3= 


137.63 


H2 


0.5000 0.0000 0.5000 



Immm 




a=2.142 


6=1.203 c= 


=2.395 


HI 


0.5000 


0.0000 


0.6766 


3 




a 


=/3= 7 =90.00 




H2 


0.0000 


0.0000 


0.5000 


C2/m(2) 


1000 


a=3.180 


6=1.101 c-- 


=2.711 


HI 


0.0006 


0.0000 


0.3376 


3 




a=7= 


90.00 0= 


=147.51 


H2 


0.5000 


0.0000 


0.0000 



Immm 




a=1.708 6=1.101 c= 


=2.711 


HI 


0.0000 


0.0000 


0.0000 


3 




a =/3= 7 =90.00 




H2 


0.5000 


0.0000 


0.1674 


Pmmn 


1104 


a=2.349 6=1.110 c= 


=1.877 


HI 


0.8170 


0.5000 


0.6886 


6 




a=/3= 7 =90.00 




H2 


0.0000 


0.5000 


0.1087 


/3-Hg (IA/mmm) 


1012 


o=6=1.408 c= 


=0.856 


HI 


0.0000 


0.0000 


0.0000 


1 




a=/ 3= 7 =90.00 












Fddd 


533 


a=3.042 6=1.832 c= 


=1.596 


HI 


0.7500 


0.7500 


0.2500 
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a=/3=7=90.00 
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